% Algorithm 1

for k=1:100
    for i=1:100
        for j=1:100
           %S1(i,j) = (i*cos(k/100*pi)+j*sin(k/100*pi))/200;
            S1(i,j) = (i+j)/200;
            S2(i,j) = (i*cos(k/100*pi)+j*sin(k/100*pi))/200;
        end
    end

    % Do QR decomposition on S1
    [Q R] = qr(S1);
    
    % Compute QTS2 to be used in CS decomposition
    QTS2 = transpose(Q)*S2;

    % CS decomposition is equivalent to four SVDs,
    % so blockify the matrix
    QTS2_1  = QTS2(1:50  , 1:50  );
    QTS2_2  = QTS2(51:100, 1:50  );
    QTS2_3  = QTS2(1:50   ,51:100);
    QTS2_4  = QTS2(51:100 ,51:100);
    
    % Compute four SVD decompositions
    [U1 E1 V1] = svd(QTS2_1);
    [U2 E2 V2] = svd(QTS2_2);
    [U3 E3 V3] = svd(QTS2_3);
    [U4 E4 V4] = svd(QTS2_4);
   
    % According to formula, velocity matrix is
    % V2*THETA*transpose(V1)
    A = V2*E1*transpose(V1);
    
   %imagesc(S1);                        % Source domain
   %imagesc(S2);                        % Target domain
   %imagesc(V1);                        % Unitary matrix
   %imagesc(U1);                        % Unitary matrix
   %imagesc(E1);                        % Angles
   imagesc(A);                         % Velocity matrix!
    colormap(flipud(gray));
    axis equal;
    M1(k) = getframe;
end

movie(M1,1,25);
